※ 본 문서는 PC 환경에 최적화되어있습니다.

1 Introduction

도시지리학에서 해당 지역의 치안은 지역에 대한 사람들의 선호를 결정하는데 결정적인 요소로서 연구대상으로 많이 사용되고 있습니다.
도시의 특성과 성격을 분석하는데 유의미한 지표라고 할 수 있는 범죄율은, 여러 유발 요인이 존재합니다. 예를 들어, 공공재로서 치안의 공급도를 한 눈에 볼 수 있는 1인당 CCTV나 경찰력, 그리고 소득 등이 있습니다.
이에 본 조는, 이러한 요소들을 활용해 범죄율에 대한 여러 분석을 진행해보기로 했습니다.

2 Methods

2.1 라이브러러리 설치

library(car)
library(DescTools)
library(GGally)
library(gridExtra)
library(ggplot2)
library(ggformula)
library(tidyverse)
library(readxl)
library(broom)
library(sf)
library(raster)
library(rgdal)
library(tmap)
library(rgeos)
library(ape)
library(spgwr)
library(spdep)
library(DT)

2.2 자료 구득

범죄율과 1만명당 CCTV, 경찰1인당 담당인구수는 2016년 SBS데이터저널리즘 마부작침 팀에서 발표한 ’범죄여지도’자료를 활용했고, GRDP와 순이동, 주민등록인구대비 기초 수급자 비율은 통계청 국가통계포털에서 발표한 2016년 기준 자료를 활용했습니다.

grdp <- read.csv(file = "./data/grdp2016.csv")
movement <- read.csv(file = "./data/emigration2016.csv")
crimedata <- read.csv(file = "./data/crime2016.csv")
tbl <- read.csv("./data/기초생활수급자.csv")
population <- read.csv("./data/인구.csv")

2.3 변수 설정

구득한 자료를 토대로 범죄유발 요인의 대리 변수를 설정하였습니다. - 범죄율 : 인구 10000명당 5대범죄 발생 건수 - 주거불안정 : 순인구이동(전입자 + 전출자) - 경제적 불이익 : 주민등록 인구 중 기초생활보장수급자가 차지하는 비율 - 소득 : 시군구별 GRDP - 경찰력 : 경찰관 1명당 담당하는 인구수 - CCTV : 인구 10000명 당 CCTV 설치대수

2.4 자료 정제하기

2.4.1 기초수급자 자료 정제

tbl <- mutate(tbl, name = paste(tbl$통계시도명, tbl$통계시군구명))
tbl <- aggregate(수급자수 ~ name, tbl, sum)
tbl <- rename(tbl, 지역명 = name)

2.4.2 순이동 자료 정제

순이동 자료 행렬에서 항목 이름 변경하기

movement <- rename(movement, 지역명 = 행정구역.시군구.별)
movement <- rename(movement, 유입인구 = X2016)
movement <- rename(movement, 유출인구 = X2016.1)
movement <- rename(movement, 순이동 = X2016.2)

순이동자료에서 합계행, 신생자치단체, 없어진 자치단체로서 유의미한 통계자료 없는 행 제거

movement <- movement[-c(1, 2, 3, 29, 46, 55, 66, 72, 78, 84, 86, 115, 116, 117, 118, 
    119, 121, 124, 125, 126, 127, 128, 147, 151, 160, 170, 171, 178, 179, 194, 199, 
    204, 219, 243, 244, 246, 247, 249, 259, 267, 270, 271), ]

2.4.3 GRDP 자료 정제

소계 없는 행렬데이터 만들고 병합

grdpdata <- subset(grdp, !행정구역소 == "소계")
grdpdata <- mutate(grdpdata, 지역명 = paste(grdpdata$행정구역대, grdpdata$행정구역소))
grdpdata <- subset(grdpdata, select = -행정구역대)
grdpdata <- subset(grdpdata, select = -행정구역소)
grdpdata <- rename(grdpdata, GRDP = X2016)

2.4.4 범죄율 자료 정제

crimedata <- mutate(crimedata, 지역명 = paste(crimedata$시도, crimedata$시군구))
crimedata <- dplyr::select(crimedata, 4, 13:15)
crimedata <- rename(crimedata, 만명당CCTV = X1만명당.CCTV)
crimedata <- rename(crimedata, 경찰한명담당인구 = X1인당.담당인구)

2.4.5 자료 통합 후 수급률 구하기

finaldata <- list(tbl, grdpdata, crimedata, movement, population) %>% reduce(inner_join, 
    by = "지역명")

t.rate <- finaldata$수급자수/finaldata$주민등록인구수
finaldata <- finaldata %>% mutate(t.rate)
finaldata <- finaldata %>% rename(crime = 만명당.5대범죄, cctv = 만명당CCTV, police = 경찰한명담당인구, 
    migration = 순이동)
finaldata$migration <- as.numeric(finaldata$migration)

2.4.6 데이터 분포 확인

2.5 datatable

str(finaldata)

'data.frame':   188 obs. of  11 variables:
 $ 지역명        : chr  "강원도 강릉시" "강원도 고성군" "강원도 동해시" "강원도 삼척시" ...
 $ 수급자수      : int  8722 1355 3623 3343 1178 1443 11630 1219 1423 1926 ...
 $ GRDP          : int  4854518 951171 2616335 2273273 867428 1297094 7945326 1381148 1329883 1592046 ...
 $ crime         : num  126.9 67.4 100.4 120.7 54.6 ...
 $ cctv          : num  30.9 67.1 53.6 31 77 ...
 $ police        : num  543 253 407 319 238 ...
 $ 유입인구      : Factor w/ 248 levels "0","10058","10358",..: 92 161 4 236 98 144 189 131 117 164 ...
 $ 유출인구      : Factor w/ 248 levels "0","10020","10131",..: 98 163 11 244 105 146 177 157 128 189 ...
 $ migration     : num  99 166 115 3 24 228 202 114 62 128 ...
 $ 주민등록인구수: int  213846 30114 93297 69599 24010 40073 337979 32720 38718 48013 ...
 $ t.rate        : num  0.0408 0.045 0.0388 0.048 0.0491 ...

head(finaldata)

         지역명 수급자수    GRDP  crime  cctv police 유입인구 유출인구
1 강원도 강릉시     8722 4854518 126.91 30.91 542.76    30907    31404
2 강원도 고성군     1355  951171  67.41 67.08 253.06     5011     4816
3 강원도 동해시     3623 2616335 100.43 53.59 407.41    10520    11151
4 강원도 삼척시     3343 2273273 120.69 31.03 319.26     8483     9508
5 강원도 양구군     1178  867428  54.56 77.05 237.72     3242     3384
6 강원도 영월군     1443 1297094  60.14 75.11 288.29     4385     4316
  migration 주민등록인구수     t.rate
1        99         213846 0.04078636
2       166          30114 0.04499568
3       115          93297 0.03883297
4         3          69599 0.04803230
5        24          24010 0.04906289
6       228          40073 0.03600928

summary(finaldata)

    지역명             수급자수          GRDP              crime       
 Length:188         Min.   :  269   Min.   :  244797   Min.   : 34.91  
 Class :character   1st Qu.: 2401   1st Qu.: 1368332   1st Qu.: 69.37  
 Mode  :character   Median : 4868   Median : 3140862   Median : 87.47  
                    Mean   : 6507   Mean   : 5653618   Mean   : 97.25  
                    3rd Qu.: 8802   3rd Qu.: 6553801   3rd Qu.:112.06  
                    Max.   :24889   Max.   :55136745   Max.   :409.00  
      cctv            police         유입인구      유출인구     migration     
 Min.   : 10.41   Min.   :125.4   10058  :  1   10131  :  1   Min.   :  2.00  
 1st Qu.: 29.68   1st Qu.:323.4   10358  :  1   10181  :  1   1st Qu.: 59.75  
 Median : 42.98   Median :447.5   10520  :  1   10267  :  1   Median :117.50  
 Mean   : 50.78   Mean   :457.2   10559  :  1   10825  :  1   Mean   :119.73  
 3rd Qu.: 64.51   3rd Qu.:581.1   105735 :  1   108949 :  1   3rd Qu.:178.25  
 Max.   :206.70   Max.   :989.8   10612  :  1   10948  :  1   Max.   :242.00  
 주민등록인구수        t.rate       
 Min.   :  10001   Min.   :0.01085  
 1st Qu.:  52135   1st Qu.:0.03029  
 Median : 113943   Median :0.04060  
 Mean   : 193754   Mean   :0.04061  
 3rd Qu.: 311200   3rd Qu.:0.04971  
 Max.   :1194041   Max.   :0.08584  
 [ getOption("max.print") 에 도달했습니다 -- 1 행들을 생략합니다 ]

DescTools::Desc(finaldata)

------------------------------------------------------------------------------ 
Describe finaldata (data.frame):

data frame: 188 obs. of  11 variables
        188 complete cases (100.0%)

  Nr  ColName         Class      NAs  Levels                                
  1   지역명          character  .                                          
  2   수급자수        integer    .                                          
  3   GRDP            integer    .                                          
  4   crime           numeric    .                                          
  5   cctv            numeric    .                                          
  6   police          numeric    .                                          
  7   유입인구        factor     .    (248): 1-0, 2-10058, 3-10358, 4-10520,
                                      5-10559, ...                          
  8   유출인구        factor     .    (248): 1-0, 2-10020, 3-10131, 4-10181,
                                      5-10267, ...                          
  9   migration       numeric    .                                          
  10  주민등록인구수  integer    .                                          
  11  t.rate          numeric    .                                          


------------------------------------------------------------------------------ 
1 - 지역명 (character)

  length      n    NAs unique levels  dupes
     188    188      0    188    188      n
         100.0%   0.0%                     

            level  freq  perc  cumfreq  cumperc
1   강원도 강릉시     1  0.5%        1     0.5%
2   강원도 고성군     1  0.5%        2     1.1%
3   강원도 동해시     1  0.5%        3     1.6%
4   강원도 삼척시     1  0.5%        4     2.1%
5   강원도 양구군     1  0.5%        5     2.7%
6   강원도 영월군     1  0.5%        6     3.2%
7   강원도 원주시     1  0.5%        7     3.7%
8   강원도 인제군     1  0.5%        8     4.3%
9   강원도 정선군     1  0.5%        9     4.8%
10  강원도 철원군     1  0.5%       10     5.3%
11  강원도 춘천시     1  0.5%       11     5.9%
12  강원도 태백시     1  0.5%       12     6.4%
... etc.
 [list output truncated]

------------------------------------------------------------------------------ 
2 - 수급자수 (integer)

     length         n       NAs    unique        0s       mean     meanCI
        188       188         0       187         0   6'507.19   5'745.04
               100.0%      0.0%                0.0%              7'269.34
                                                                         
        .05       .10       .25    median       .75        .90        .95
   1'327.70  1'468.00  2'400.75  4'868.50  8'801.50  14'261.30  16'709.35
                                                                         
      range        sd     vcoef       mad       IQR       skew       kurt
  24'620.00  5'297.26      0.81  4'241.72  6'400.75       1.26       1.12
                                                                         
lowest : 269, 898, 1'040, 1'047, 1'071
highest: 21'071, 22'680, 23'335, 23'783, 24'889

------------------------------------------------------------------------------ 
3 - GRDP (integer)

      length             n           NAs        unique            0s
         188           188             0           = n             0
                    100.0%          0.0%                        0.0%
                                                                    
         .05           .10           .25        median           .75
  647'101.45    794'546.30  1'368'331.75  3'140'862.00  6'553'801.00
                                                                    
       range            sd         vcoef           mad           IQR
    5.49e+07  7'666'003.83          1.36  3'232'469.04  5'185'469.25
                                                                    
          mean        meanCI
  5'653'617.78  4'550'662.29
                6'756'573.26
                            
           .90           .95
      1.09e+07      2.22e+07
                            
          skew          kurt
          3.46         15.28
                            
lowest : 244'797, 385'039, 459'213, 476'103, 548'014
highest: 3e+07, 3e+07, 3e+07, 5e+07, 6e+07

------------------------------------------------------------------------------ 
4 - crime (numeric)

    length        n      NAs   unique        0s      mean    meanCI
       188      188        0      186         0   97.2522   90.2469
             100.0%     0.0%               0.0%            104.2576
                                                                   
       .05      .10      .25   median       .75       .90       .95
   51.7365  54.6300  69.3675  87.4750  112.0625  132.7120  160.3975
                                                                   
     range       sd    vcoef      mad       IQR      skew      kurt
  374.0900  48.6902   0.5007  30.4081   42.6950    3.3009   15.7394
                                                                   
lowest : 34.91, 37.06, 38.76, 43.71, 44.14
highest: 217.77, 296.59, 306.73, 377.73, 409.0

------------------------------------------------------------------------------ 
5 - cctv (numeric)

    length        n      NAs   unique       0s     mean    meanCI
       188      188        0      = n        0  50.7777   46.3944
             100.0%     0.0%              0.0%            55.1610
                                                                 
       .05      .10      .25   median      .75      .90       .95
   16.0365  20.1300  29.6775  42.9800  64.5125  84.7670  102.5735
                                                                 
     range       sd    vcoef      mad      IQR     skew      kurt
  196.2900  30.4656   0.6000  23.4547  34.8350   1.7483    4.7513
                                                                 
lowest : 10.41, 10.89, 11.81, 12.13, 12.37
highest: 131.67, 133.12, 155.14, 177.84, 206.7

------------------------------------------------------------------------------ 
6 - police (numeric)

    length         n       NAs    unique        0s      mean    meanCI
       188       188         0       = n         0  457.2343  433.0252
              100.0%      0.0%                0.0%            481.4434
                                                                      
       .05       .10       .25    median       .75       .90       .95
  230.8050  253.5850  323.3550  447.5000  581.0625  672.9870  739.6410
                                                                      
     range        sd     vcoef       mad       IQR      skew      kurt
  864.4000  168.2636    0.3680  190.0990  257.7075    0.4412   -0.3280
                                                                      
lowest : 125.37, 136.37, 144.9, 172.43, 202.09
highest: 811.2, 827.04, 870.98, 883.27, 989.77

------------------------------------------------------------------------------ 
7 - 유입인구 (factor)

  length      n    NAs unique levels  dupes
     188    188      0    188    248      n
         100.0%   0.0%                     

     level  freq  perc  cumfreq  cumperc
1    10058     1  0.5%        1     0.5%
2    10358     1  0.5%        2     1.1%
3    10520     1  0.5%        3     1.6%
4    10559     1  0.5%        4     2.1%
5   105735     1  0.5%        5     2.7%
6    10612     1  0.5%        6     3.2%
7    10957     1  0.5%        7     3.7%
8    11446     1  0.5%        8     4.3%
9    11899     1  0.5%        9     4.8%
10   12210     1  0.5%       10     5.3%
11   12245     1  0.5%       11     5.9%
12   12402     1  0.5%       12     6.4%
... etc.
 [list output truncated]

------------------------------------------------------------------------------ 
8 - 유출인구 (factor)

  length      n    NAs unique levels  dupes
     188    188      0    188    248      n
         100.0%   0.0%                     

     level  freq  perc  cumfreq  cumperc
1    10131     1  0.5%        1     0.5%
2    10181     1  0.5%        2     1.1%
3    10267     1  0.5%        3     1.6%
4    10825     1  0.5%        4     2.1%
5   108949     1  0.5%        5     2.7%
6    10948     1  0.5%        6     3.2%
7    11151     1  0.5%        7     3.7%
8    11154     1  0.5%        8     4.3%
9   112735     1  0.5%        9     4.8%
10   11697     1  0.5%       10     5.3%
11   11819     1  0.5%       11     5.9%
12   12112     1  0.5%       12     6.4%
... etc.
 [list output truncated]

------------------------------------------------------------------------------ 
9 - migration (numeric)

  length       n    NAs  unique      0s    mean  meanCI
     188     188      0     186       0  119.73  109.81
          100.0%   0.0%            0.0%          129.66
                                                       
     .05     .10    .25  median     .75     .90     .95
   15.35   25.00  59.75  117.50  178.25  213.30  227.65
                                                       
   range      sd  vcoef     mad     IQR    skew    kurt
  240.00   68.97   0.58   88.21  118.50    0.02   -1.21
                                                       
lowest : 2.0, 3.0, 4.0, 6.0, 7.0
highest: 238.0, 239.0, 240.0, 241.0, 242.0

------------------------------------------------------------------------------ 
10 - 주민등록인구수 (integer)

        length           n        NAs      unique          0s        mean
           188         188          0         = n           0  193'754.14
                    100.0%       0.0%                    0.0%            
                                                                         
           .05         .10        .25      median         .75         .90
     27'656.65   33'303.20  52'134.75  113'942.50  311'200.00  440'587.90
                                                                         
         range          sd      vcoef         mad         IQR        skew
  1'184'040.00  178'154.75       0.92  119'788.89  259'065.25        1.56
                                                                         
      meanCI
  168'121.91
  219'386.37
            
         .95
  501'469.75
            
        kurt
        4.13
            
lowest : 10'001, 17'713, 23'628, 24'010, 24'171
highest: 591'891, 595'485, 657'831, 662'154, 1'194'041

------------------------------------------------------------------------------ 
11 - t.rate (numeric)

      length           n         NAs      unique          0s        mean
         188         188           0         = n           0  0.04061096
                  100.0%        0.0%                    0.0%            
                                                                        
         .05         .10         .25      median         .75         .90
  0.01647370  0.02050721  0.03029066  0.04060200  0.04971344  0.05797234
                                                                        
       range          sd       vcoef         mad         IQR        skew
  0.07498548  0.01476651  0.36360904  0.01407992  0.01942278  0.22213282
                                                                        
       meanCI
   0.03848641
   0.04273551
             
          .95
   0.06553814
             
         kurt
  -0.06175245
             
lowest : 0.01085216, 0.01318128, 0.01334563, 0.01342153, 0.01408009
highest: 0.07179431, 0.07321147, 0.07770532, 0.08268207, 0.08583764

DT::datatable(finaldata, rownames = FALSE, caption = "show finaldata")

2.6 대수변환

2.6.1 지역별 5대 범죄율

# crime data
par(mfrow = c(2, 2))
hist(finaldata$crime)
qqnorm(finaldata$crime)
qqline(finaldata$crime)
log.crime <- log(finaldata$crime)
hist(log.crime)
qqnorm(log.crime)
qqline(log.crime)

finaldata <- finaldata %>% mutate(log.crime)

2.6.2 지역별 GRDP

# GRDP
par(mfrow = c(2, 2))
hist(finaldata$GRDP)
qqnorm(finaldata$GRDP)
qqline(finaldata$GRDP)
log.GRDP <- log(finaldata$GRDP)
hist(log.GRDP)
qqnorm(log.GRDP)
qqline(log.GRDP)

finaldata <- finaldata %>% mutate(log.GRDP)

2.6.3 지역별 1만명당 CCTV 대수

par(mfrow = c(2, 2))
hist(finaldata$cctv)
qqnorm(finaldata$cctv)
qqline(finaldata$cctv)
log.cctv <- log(finaldata$cctv)
hist(log.cctv)
qqnorm(log.cctv)
qqline(log.cctv)

finaldata <- finaldata %>% mutate(log.cctv)

2.6.4 Summary

DT::datatable(finaldata, rownames = FALSE, caption = "FINAL DATA TABLE")

p1 <- ggplot(data = finaldata, aes(y = crime, x = GRDP)) + geom_point(alpha = 0.2)
p2 <- ggplot(data = finaldata, aes(y = log.crime, x = log.GRDP)) + geom_point(alpha = 0.2)
grid.arrange(p1, p2, ncol = 2)

로그를 씌우기 전과 로그를 씌운 후 GRDP 변수와 종속변수인 범죄율의 관계가 어떻게 나타나는지를 살펴봤습니다.
눈으로 보기에도 두 변수 간의 관계가 양의 상관관계를 보임을 알 수 있지만, 로그 값을 적용함으로써 GRDP 변수가 미치는 영향을 조금 더 선명하게 확인할 수 있습니다.

ggpairs(data = finaldata, columns = c("crime", "GRDP", "cctv", "police", "t.rate", 
    "migration"))

ggpairs(data = finaldata, columns = c("log.crime", "log.GRDP", "log.cctv", "police", 
    "t.rate", "migration"))

3 Analysis & Results

3.1 단순 선형 회귀분석

3.1.1 상관관계 분석 및 시각화

cor.test(finaldata$log.GRDP, finaldata$log.crime)


    Pearson's product-moment correlation

data:  finaldata$log.GRDP and finaldata$log.crime
t = 10.624, df = 186, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5168939 0.6964074
sample estimates:
      cor 
0.6145446

상관관계를 알아본 결과 log.GRDP와 log.crime은0.6145446으로 강한 양의 상관관계를 보였습니다.

ggplot(finaldata, aes(x = log.GRDP, y = log.crime)) + geom_point() + geom_smooth(method = "lm") + 
    theme_classic()

p1 <- ggplot(finaldata, aes(x = log.GRDP, y = log.crime)) + geom_point() + xlab("GRDP (log scale)") + 
    ylab("Crime rate(log scale)") + geom_smooth(method = lm, color = "blue", fill = "#69b3a2", 
    se = TRUE) + theme_classic()
p1

그래프를 통해서도 log.crime과 log.GRDP가 선형관계를 가짐을 알 수 있습니다. 이에 이 두 변수를 가지고 회귀분석을 진행하였습니다.

3.1.2 회귀 분석

mod1 <- lm(log.crime ~ log.GRDP, data = finaldata)
summary(mod1)


Call:
lm(formula = log.crime ~ log.GRDP, data = finaldata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57378 -0.18341 -0.03834  0.14911  1.56077 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.14477    0.31599   3.623 0.000376 ***
log.GRDP     0.22373    0.02106  10.624  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3068 on 186 degrees of freedom
Multiple R-squared:  0.3777,    Adjusted R-squared:  0.3743 
F-statistic: 112.9 on 1 and 186 DF,  p-value: < 2.2e-16

anova(mod1)

Analysis of Variance Table

Response: log.crime
           Df Sum Sq Mean Sq F value    Pr(>F)    
log.GRDP    1 10.626 10.6258  112.87 < 2.2e-16 ***
Residuals 186 17.510  0.0941                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

분석 결과를 자세히 보겠습니다.
Residuals는 \(y_i - \hat{y_i}\) 값을 나타내고
Coefficients는 추정된 회귀식을 나타내며,
Estimate는 추정된 회귀식 절편과 기울기를 보여줍니다.

도출된 회귀식은 \[ \hat{y_i}= 1.14477 + 0.22373x_i \]

즉, GRDP가 1% 증가하면 범죄율은는 약0.2%만큼 증가합니다.

\(p-value < 2.2e-16\)이므로 귀무가설(\(H_0\):회귀식은 유의하지 않다)을 기각합니다. \(R^2\)는 0.3743으로 약 37%의 설명력을 가진다고 볼 수 있습니다.

tidy(mod1)

# A tibble: 2 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    1.14     0.316       3.62 3.76e- 4
2 log.GRDP       0.224    0.0211     10.6  6.59e-21

3.1.3 분석 검정

그렇다면, 잔차 분석을 통해 회귀식이 적합한지 확인 해보겠습니다.

par(mfrow = c(2, 2))
plot(mod1)

durbinWatsonTest(mod1)

 lag Autocorrelation D-W Statistic p-value
   1      0.06372597      1.868454    0.36
 Alternative hypothesis: rho != 0

plot(scale(mod1$fitted), scale(mod1$res))
abline(h = 0)

3.2 다중회귀분석

확장판으로 고고! 다중 회귀 분석은 예측자가 2개 이상인 경우에 진행합니다.

\[ \hat{y_i} = b_0 + b_1x_{1i} + b_2x_{2i} \]

회귀식이 위와 같을 때,
\(x_{1i}\)가 1 단위 증가하면 \(\hat{y_i}\)는 \(b_1\)만큼 변하며(단,\(x_{2i}\)는 고정)
\(x_{2i}\)가 1 단위 증가하면 \(\hat{y_i}\)는 \(b_2\)만큼 변한다고 해석합니다(단,\(x_{1i}\)는 고정)

3.2.1 회귀 분석

mod4 <- lm(log.crime ~ log.GRDP + log.cctv + police + migration + t.rate, data = finaldata)
summary(mod4)


Call:
lm(formula = log.crime ~ log.GRDP + log.cctv + police + migration + 
    t.rate, data = finaldata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.62652 -0.17008 -0.02112  0.15570  1.18872 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.6951759  0.5248531  -1.325 0.186993    
log.GRDP     0.3553564  0.0278827  12.745  < 2e-16 ***
log.cctv     0.0012125  0.0416701   0.029 0.976819    
police      -0.0007838  0.0001729  -4.534 1.05e-05 ***
migration   -0.0003017  0.0003040  -0.993 0.322174    
t.rate       6.3944769  1.8456417   3.465 0.000662 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.27 on 182 degrees of freedom
Multiple R-squared:  0.5284,    Adjusted R-squared:  0.5155 
F-statistic: 40.79 on 5 and 182 DF,  p-value: < 2.2e-16

anova(mod4)

Analysis of Variance Table

Response: log.crime
           Df  Sum Sq Mean Sq  F value    Pr(>F)    
log.GRDP    1 10.6258 10.6258 145.7524 < 2.2e-16 ***
log.cctv    1  0.0088  0.0088   0.1202 0.7292397    
police      1  3.2488  3.2488  44.5628 2.881e-10 ***
migration   1  0.1087  0.1087   1.4914 0.2235853    
t.rate      1  0.8751  0.8751  12.0037 0.0006621 ***
Residuals 182 13.2684  0.0729                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(mod4)

                    2.5 %        97.5 %
(Intercept) -1.7307550825  0.3404033557
log.GRDP     0.3003415221  0.4103712293
log.cctv    -0.0810061879  0.0834311012
police      -0.0011248429 -0.0004426676
migration   -0.0009014686  0.0002979940
t.rate       2.7528705495 10.0360832058

tidy(mod4)

# A tibble: 6 x 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) -0.695     0.525      -1.32   1.87e- 1
2 log.GRDP     0.355     0.0279     12.7    5.28e-27
3 log.cctv     0.00121   0.0417      0.0291 9.77e- 1
4 police      -0.000784  0.000173   -4.53   1.05e- 5
5 migration   -0.000302  0.000304   -0.993  3.22e- 1
6 t.rate       6.39      1.85        3.46   6.62e- 4

다중회귀분석 시행 결과
\(F-statistic)\)의 \(p-value\)가 0에 가까운 값이므로 귀무가설(\(H_0\):회귀식은 유의하지 않다)을 기각합니다. \(R^2\)는 0.5284로 약 52%의 설명력을 가진다고 볼 수 있습니다.

하지만, 일부 변수들의 경우 p값이 0.05를 초과해 아래의 변수 선택 과정을 통해 유의미한 변수들을 선택해보겠습니다.

3.2.2 변수 선택

3.2.2.1 전진선택법

# forward selection
summary(step(mod4, direction = "forward"))

Start:  AIC=-486.4
log.crime ~ log.GRDP + log.cctv + police + migration + t.rate


Call:
lm(formula = log.crime ~ log.GRDP + log.cctv + police + migration + 
    t.rate, data = finaldata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.62652 -0.17008 -0.02112  0.15570  1.18872 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.6951759  0.5248531  -1.325 0.186993    
log.GRDP     0.3553564  0.0278827  12.745  < 2e-16 ***
log.cctv     0.0012125  0.0416701   0.029 0.976819    
police      -0.0007838  0.0001729  -4.534 1.05e-05 ***
migration   -0.0003017  0.0003040  -0.993 0.322174    
t.rate       6.3944769  1.8456417   3.465 0.000662 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.27 on 182 degrees of freedom
Multiple R-squared:  0.5284,    Adjusted R-squared:  0.5155 
F-statistic: 40.79 on 5 and 182 DF,  p-value: < 2.2e-16

독립변수 중 종속변수에 가장 큰 영향을 주는 변수부터 모형에 포함하는 전진선택법을 적용한 결과입니다.

3.2.2.2 후진제거법

# backward elimination
backward <- step(mod4, direction = "backward")

Start:  AIC=-486.4
log.crime ~ log.GRDP + log.cctv + police + migration + t.rate

            Df Sum of Sq    RSS     AIC
- log.cctv   1    0.0001 13.268 -488.40
- migration  1    0.0718 13.340 -487.38
<none>                   13.268 -486.40
- t.rate     1    0.8751 14.143 -476.39
- police     1    1.4985 14.767 -468.28
- log.GRDP   1   11.8415 25.110 -368.48

Step:  AIC=-488.4
log.crime ~ log.GRDP + police + migration + t.rate

            Df Sum of Sq    RSS     AIC
- migration  1    0.0736 13.342 -489.36
<none>                   13.268 -488.40
- t.rate     1    0.9355 14.204 -477.59
- police     1    1.6996 14.968 -467.74
- log.GRDP   1   12.8049 26.073 -363.40

Step:  AIC=-489.36
log.crime ~ log.GRDP + police + t.rate

           Df Sum of Sq    RSS     AIC
<none>                  13.342 -489.36
- t.rate    1    1.0141 14.356 -477.59
- police    1    1.9397 15.282 -465.84
- log.GRDP  1   14.0654 27.407 -356.02

summary(backward)


Call:
lm(formula = log.crime ~ log.GRDP + police + t.rate, data = finaldata)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6177 -0.1647 -0.0211  0.1552  1.1817 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.8107103  0.4023128  -2.015 0.045347 *  
log.GRDP     0.3615319  0.0259581  13.928  < 2e-16 ***
police      -0.0008201  0.0001586  -5.172 6.01e-07 ***
t.rate       6.5950486  1.7635095   3.740 0.000246 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2693 on 184 degrees of freedom
Multiple R-squared:  0.5258,    Adjusted R-squared:  0.5181 
F-statistic: 68.01 on 3 and 184 DF,  p-value: < 2.2e-16

위와 반대로, 독립변수를 모두 포함한 모형에서 가장 영향이 적은 변수부터 제거하는 후진 제거법을 적용한 결과입니다.

3.2.2.3 단계적 회귀

# stepwise regression
summary(step(mod4, direction = "both"))

Start:  AIC=-486.4
log.crime ~ log.GRDP + log.cctv + police + migration + t.rate

            Df Sum of Sq    RSS     AIC
- log.cctv   1    0.0001 13.268 -488.40
- migration  1    0.0718 13.340 -487.38
<none>                   13.268 -486.40
- t.rate     1    0.8751 14.143 -476.39
- police     1    1.4985 14.767 -468.28
- log.GRDP   1   11.8415 25.110 -368.48

Step:  AIC=-488.4
log.crime ~ log.GRDP + police + migration + t.rate

            Df Sum of Sq    RSS     AIC
- migration  1    0.0736 13.342 -489.36
<none>                   13.268 -488.40
+ log.cctv   1    0.0001 13.268 -486.40
- t.rate     1    0.9355 14.204 -477.59
- police     1    1.6996 14.968 -467.74
- log.GRDP   1   12.8049 26.073 -363.40

Step:  AIC=-489.36
log.crime ~ log.GRDP + police + t.rate

            Df Sum of Sq    RSS     AIC
<none>                   13.342 -489.36
+ migration  1    0.0736 13.268 -488.40
+ log.cctv   1    0.0018 13.340 -487.38
- t.rate     1    1.0141 14.356 -477.59
- police     1    1.9397 15.282 -465.84
- log.GRDP   1   14.0654 27.407 -356.02


Call:
lm(formula = log.crime ~ log.GRDP + police + t.rate, data = finaldata)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6177 -0.1647 -0.0211  0.1552  1.1817 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.8107103  0.4023128  -2.015 0.045347 *  
log.GRDP     0.3615319  0.0259581  13.928  < 2e-16 ***
police      -0.0008201  0.0001586  -5.172 6.01e-07 ***
t.rate       6.5950486  1.7635095   3.740 0.000246 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2693 on 184 degrees of freedom
Multiple R-squared:  0.5258,    Adjusted R-squared:  0.5181 
F-statistic: 68.01 on 3 and 184 DF,  p-value: < 2.2e-16

마지막으로, 전진선택법에 의해 변수를 추가하고, 변수 추가시 기존의 변수 중요도가 정해진 유의수준에 포함되지 않으면 앞에서 들어간 변수도 다시 제거하는 단계적 회귀방법을 적용하였습니다.

3.2.3 fitting

결정계수 값이 조금 더 큰 후진제거법을 적용한 모델을 선택하였습니다.

fit <- backward
summary(fit)


Call:
lm(formula = log.crime ~ log.GRDP + police + t.rate, data = finaldata)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6177 -0.1647 -0.0211  0.1552  1.1817 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.8107103  0.4023128  -2.015 0.045347 *  
log.GRDP     0.3615319  0.0259581  13.928  < 2e-16 ***
police      -0.0008201  0.0001586  -5.172 6.01e-07 ***
t.rate       6.5950486  1.7635095   3.740 0.000246 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2693 on 184 degrees of freedom
Multiple R-squared:  0.5258,    Adjusted R-squared:  0.5181 
F-statistic: 68.01 on 3 and 184 DF,  p-value: < 2.2e-16

3.3 잔차 분석

3.3.1 Residusal

par(mfrow = c(2, 2))
plot(fit)

3.3.2 정규성 검정

qqnorm(fit$residuals)
qqline(fit$residuals)

3.3.3 독립성 검정

durbinWatsonTest(fit)

 lag Autocorrelation D-W Statistic p-value
   1      0.05403652      1.886683   0.404
 Alternative hypothesis: rho != 0

독립성을 Durbin-Watson테스트를 통해 검정해본 결과, p-value 값이 0.448로 자기상관이 없다고 판단을 내릴 수 있습니다.

3.3.4 선형성 검정

crPlots(fit)

선형성을 검정한 결과, 변수들이 모두 선형 관계를 보이므로 선형성가정에 부합합니다.

3.3.5 등분산성 검정

plot(fit, which = 1)

plot(fit, which = 3)

3.3.6 다중공선성 검정

독립변수 간의 상관관계가 나타나는 다중공선성이 존재한다면 유의미한 회귀분석이라고 보기 어렵습니다.
만약 \(x_{1i}\)과 \(x_{2i}\)의 상관성이 높은 경우라면, \(x_{1i}\)이 증가하였을 때,\(x_{2i}\)는 함께 증가하는 경우가 발생해버리는데, 이를 다중공선성(Multicollinearity)이라고 합니다.

다중공선성은 분산팽창지수인(VIF)를 통해 검정할 수 있습니다.

# Multicollinearity
vif(mod4)

 log.GRDP  log.cctv    police migration    t.rate 
 2.263832  1.458784  2.170284  1.127152  1.905225

sqrt(vif(mod4)) > 2

 log.GRDP  log.cctv    police migration    t.rate 
    FALSE     FALSE     FALSE     FALSE     FALSE

보통, VIF가 2 이상일 때 다중공선성의 존재를 의심해 볼 수 있는데, 검정 결과 2가 넘는 값이 나타나지 않기 때문에, 다중공선성은 존재하지 않습니다.

3.4 Influence Statistics

3.4.1 Residuals graph

qqPlot(fit, labels = row.names(지역명), id.method = "identify", simulate = TRUE)

[1]  86 103

residplot <- function(finaldata, nbreaks = 10) {
    z <- rstudent(fit)
    hist(z, breaks = nbreaks, freq = FALSE, xlab = "Studentized Residual", main = "Distribution of Errors")
    rug(jitter(z), col = "brown")
    curve(dnorm(x, mean = mean(z), sd = sd(z)), add = TRUE, col = "blue", lwd = 2)
    lines(density(z)$x, density(z)$y, col = "red", lwd = 2, lty = 2)
    legend("topright", legend = c("Normal Curve", "Kernel Density Curve"), lty = 1:2, 
        col = c("blue", "red"), cex = 0.7)
}
residplot(fit)

3.4.2 Hat Matrix

inf.fit <- influence.measures(fit)
print(round(inf.fit$infmat, 4), 3)

     dfb.1_ dfb.l.GR dfb.polc dfb.t.rt   dffit cov.r cook.d    hat
1   -0.0267   0.0193   0.0322   0.0320  0.0874 1.008 0.0019 0.0078
2   -0.0124   0.0087   0.0101   0.0088 -0.0219 1.039 0.0001 0.0171
3    0.0131  -0.0058  -0.0144  -0.0160  0.0470 1.021 0.0006 0.0066
4   -0.0107   0.0220  -0.0508   0.0114  0.0937 1.011 0.0022 0.0094
5   -0.0638   0.0459   0.0572   0.0318 -0.1282 1.018 0.0041 0.0168
6   -0.0588   0.0373   0.0488   0.0646 -0.0972 1.032 0.0024 0.0198
7   -0.0319   0.0246   0.0416   0.0204  0.0924 1.018 0.0021 0.0118
8   -0.0736   0.0409   0.0828   0.0876 -0.1416 1.020 0.0050 0.0194
9    0.1407  -0.0735  -0.1752  -0.1769  0.2803 0.971 0.0194 0.0226
 [ getOption("max.print") 에 도달했습니다 -- 179 행들을 생략합니다 ]

hat.plot <- function(fit) {
    p <- length(coefficients(fit))
    n <- length(fitted(fit))
    
    y = hatvalues(fit)
    
    name = attr(hatvalues(fit), "name")
    df = data.frame(x = 1:length(y), y = as.numeric(y), name = name)
    
    require(ggplot2)
    require(ggiraph)
    require(moonBook2)
    p1 <- ggplot(df, aes(x = x, y = y, tooltip = name, data_id = x)) + geom_point_interactive()
    yintercept2 = 2 * p/n
    p1 <- p1 + geom_hline(aes(yintercept = yintercept2), col = "red", lty = "dashed")
    yintercept3 = 3 * p/n
    p1 <- p1 + geom_hline(aes(yintercept = yintercept3), col = "red", lty = "dashed")
    ggiraph(code = print(p1))
    
}
hat.plot(fit)

3.4.3 Cook’s distance

plot(fit, which = 4)

cooksd <- cooks.distance(fit)
plot(cooksd, pch = "*", cex = 2, main = "Influential Obs by Cooks Distance")
abline(h = 4 * mean(cooksd, na.rm = T), col = "red")
text(x = 1:length(cooksd) + 1, y = cooksd, labels = ifelse(cooksd > 4 * mean(cooksd, 
    na.rm = T), names(cooksd), ""), col = "red")

Cook’s distance를 통해 아웃라이어의 존재에 대해서 확인해본 결과, 기준선의 위쪽에 일부 수치들이 존재하지만, 거의 대부분의 수치가 아래쪽에 위치함을 확인할 수 있습니다.

3.4.4 influence plots

influenceIndexPlot(fit, vars = c("Studentized", "hat", "Cook"), main = "Influence Index plot")

위의 시각화를 통해 어떤 데이터가 통계에 가장 큰 영향을 주는 아웃라이어인지를 눈으로 교차검증할 수 있습니다.

3.4.5 outlerTest

car::outlierTest(fit)

    rstudent unadjusted p-value Bonferroni p
103 4.731996         4.4346e-06    0.0008337

finaldata["103", -c(6, 7, 8, 10)]

             지역명 수급자수    GRDP crime  cctv migration     t.rate log.crime
103 부산광역시 중구     2850 2640768   409 84.28        97 0.06304194  6.013715
    log.GRDP log.cctv
103 14.78658 4.434145

fitted(fit)["103"]

     103 
4.832047

회귀식에 103번 데이터를 적용해본 결과, 실제값에 비해 적은 수치로 추정됨을 알 수 있었습니다.

3.4.6 presumption

outlier test결과 103번 데이터는 부산광역시 중구였습니다. 또한 플롯 상에서 두번째로 눈에 띄는 데이터는 86번 데이터로, 대구광역시 중구였습니다.
뿐만 아니라, 만 명당 5대 범죄율 데이터에서 서울 특별시 중구와 종로구가 높은 순위를 차지한 것을 확인할 수 있었습니다. 즉, ‘구도심’ 지역이 높은 범죄율을 보임을 알 수 있습니다.
이러한 점으로 미루어보아, 본 조의 연구에서 적용한 독립 변수 이외에 해당 지역의 범죄율에 영향을 미치는 또 다른 독립변수가 존재 할 것으로 추정하였습니다.
가령, 구도심 지역은 거주인구 대비 주간인구가 많은 지역이므로 외부인의 출입이 잦다는 특징이 있기 때문에 외부인의 유입 정도가 범죄율에 큰 영향을 미치는 변수일 가능성을 생각해 볼 수 있습니다.

3.5 Refining Model

아웃라이어였던 103번 데이터를 제외하고 분석을 진행해 봤습니다.

re = finaldata[-c(103), ]
remod4 = lm(log.crime ~ log.GRDP + police + t.rate, data = re)
summary(remod4)


Call:
lm(formula = log.crime ~ log.GRDP + police + t.rate, data = re)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.60126 -0.15690 -0.01363  0.16049  0.89696 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.6187409  0.3829412  -1.616 0.107869    
log.GRDP     0.3472000  0.0247551  14.025  < 2e-16 ***
police      -0.0007242  0.0001514  -4.782 3.55e-06 ***
t.rate       5.9098790  1.6754178   3.527 0.000531 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2549 on 183 degrees of freedom
Multiple R-squared:  0.5395,    Adjusted R-squared:  0.5319 
F-statistic: 71.45 on 3 and 183 DF,  p-value: < 2.2e-16

anova(remod4)

Analysis of Variance Table

Response: log.crime
           Df  Sum Sq Mean Sq F value    Pr(>F)    
log.GRDP    1 10.7517 10.7517 165.516 < 2.2e-16 ***
police      1  2.3646  2.3646  36.401  8.71e-09 ***
t.rate      1  0.8083  0.8083  12.443 0.0005305 ***
Residuals 183 11.8875  0.0650                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidy(remod4)

# A tibble: 4 x 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) -0.619     0.383        -1.62 1.08e- 1
2 log.GRDP     0.347     0.0248       14.0  8.04e-31
3 police      -0.000724  0.000151     -4.78 3.55e- 6
4 t.rate       5.91      1.68          3.53 5.31e- 4

이상치로 의심되는 데이터를 제외하고 분석을 진행해본 결과, p값은 이전과 같이 기각역보다 낮게 나타나 귀무가설이 기각되고, 설명계수는 소폭의 상승을 보였습니다.
변수 간 상관관계는 아웃라이어가 제거돼, 일부 강화된 것을 볼 수 있습니다.

3.6 Prediction

회귀식이 예측을 잘 하는지 확인해보겠습니다.

newdata <- re %>% dplyr::select(log.crime, log.GRDP, police, t.rate)
head(predict(remod4, newdata = newdata))

       1        2        3        4        5        6 
4.574512 4.243270 4.446375 4.515782 4.246418 4.272345

pre <- as.data.frame(predict(remod4, newdata = newdata, interval = "confidence"))

head(predict(remod4, newdata = newdata, interval = "confidence"))

       fit      lwr      upr
1 4.574512 4.530049 4.618976
2 4.243270 4.177295 4.309246
3 4.446375 4.405558 4.487192
4 4.515782 4.466199 4.565364
5 4.246418 4.180945 4.311892
6 4.272345 4.201484 4.343205

ppre <- cbind(pre, re$log.crime)
head(ppre)

       fit      lwr      upr re$log.crime
1 4.574512 4.530049 4.618976     4.843478
2 4.243270 4.177295 4.309246     4.210793
3 4.446375 4.405558 4.487192     4.609461
4 4.515782 4.466199 4.565364     4.793225
5 4.246418 4.180945 4.311892     3.999301
6 4.272345 4.201484 4.343205     4.096675

3.7 결과

선형 회귀 분석 결과, GRDP, 기초수급자 비율, 경찰력 변수가 범죄율에 영향을 미치는 것으로 나타났습니다.
가장 뚜렷한 상관관계를 보이는 것은 GRDP와 범죄율이었는데, 이는 ‘부유한 지역일수록 범죄가 적게 일어날 것’이라는 사회적 통념과 반하는 결과로, 오히려 GRDP와 범죄율 사이에는 양의 상관관계가 나타나며, 부유한 지역이 범죄로부터의 안전지대가 아니라는 것을 확인할 수 있습니다. GRDP가 높은 지역의 경우, 기업을 비롯한 각종 사무실, 고위기관들이 밀집되어 유동인구가 많고, 유흥업소 수 또한 많아 범죄율이 높아지는 결과가 도출된 것으로 추론해볼 수 있습니다.
두번째로 영향을 미치는 변수는 기초생활 수급 비율로, 이 역시 범죄율과 양의 관계를 보였는데, 그 지역의 전반적인 소득 수준을 나타내는 지표인 GRDP와 인구 대비 기초 생활 수급자 비율이 동시에 양의 관계를 보이고 있으므로 소득 격차가 심한 지역일수록 범죄율이 높게 나타나는 것으로 해석 할 수 있습니다.
마지막으로 경찰 1명당 인구 수와 범죄율은 음의 관계를 보이고 있는데, 이는 범죄율이 높은 지역에 경찰인력을 그만큼 많이 배치된 경향이 반영된 것으로 추정해보았습니다.

4 Geographically Weighted Regression(GWR)

종래의 일반최소제곱법(Ordinary Least Squares, OLS) 및 공간회귀모형은 종속변수에 대한 독립변수들의 영향력이 모든 공간 영역 안에서 모두 동일하다고 가정합니다. 즉, 기존의 공간회귀모형에서는 추정된 회귀계수가 모든 연구 대상 지역에서 동일하게 적용된다는 의미입니다. 이는 변수들 간의 독립성과 오차의 동분산성(homoscedasticity)을 가정합니다(Berry and Feldman, 1985; 김두섭·강남준, 2008). 그러나 이러한 공간적 정상성(spatial stationarity)은 국지적 지역 및 공간에서 일어나는 개별적 특성에 대한 고려가 불가능하며 지역별 다양한 요인으로 인해 영향을 받는 인구이동에 대한 분석을 진행하기에는 한계가 있습니다. 따라서 기존의 OLS회귀모형으로 공간적 비정상성 또는 의존성이란 특성들이 존재하는 현상을 분석하게 되면 모수 추정값의 정확성 및 효율성이 낮아지게 됩니다(조동기, 2009).

지리적가중회귀모형(GWR)은 이러한 공간적 비정상성을 고려하기 위해 고안된 공간분석 기법으로서, 기존 모형들의 장점들이 공유됨과 함께 지역별로 각기 다른 회귀계수의 추정이 가능하다는 장점이 있습니다. 즉, OLS모형으로 파악할 수 없었던 복잡한 공간적 이질성을 나타낼 수 있다는 것입니다. 그 뿐만 아니라 지리적가중회귀분석의 변수별 결과를 GIS 프로그램을 통해 시각적으로 제시할 수 있어 각 변수별 공간적 특성 및 상호작용을 이해할 수 있습니다(Fotheringham et al., 2002).

변수들에 내재한 공간적 비정상성을 고려하기 위해 지리적가중회귀모형은 회귀계수를 하나의 상수가 아니라, 공간적 위치에 따라 달라지도록 모형화하는 것입니다. 따라서 지리적 위치 i에 따른 변수 k에 대해서 종속변수 Y에 영향을 미치는 개별적인 회귀계수 \(\beta_i\)를 추정하는 함수식은 식 (1)와 같이 나타낼 수 있습니다.

\(Y_i = \beta_{0i} + \sum_{k=1}^m \beta_{ki} X_{ki} + \varepsilon_i\) …………………………………………………(1)

회귀계수 \(\beta_i\)는 ‘지리학 제1법칙’에 의해서 위치에 따라 가중치가 다르게 부과되며, 지리적 위치 i와 주변 지역 간 거리로 측정된 공간행렬 \(W_i\)를 바탕으로 가중회귀초소제곱(weighted least square)의 수식 (2)에 의해 회귀계수가 결정됩니다.

\(\hat{\beta}_i = (\hat{\beta}_{i0},\hat{\beta}_{i1},\hat{\beta}_{i2}, ..., \hat{\beta}_{im})= (X^TW_iX)^{-1}(X^TW_iY)\) ………………(2)

GWR은 행렬연산에 기반하기 때문에 본 연구에서는 결측치가 존재하지 않는 서울시로 한정하여 진행하였습니다.

4.1 셰이프 만들기

위의 분석 데이터를 공간상에 표현하기 위해 ETL에서 구득한 셰이프 파일을 불러왔습니다.

`2017shp` <- readOGR("./data/shapefiles", layer = "kr_si_gun-gu", encoding = "UTF-8")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Zulu-06\data\shapefiles", layer: "kr_si_gun-gu"
with 250 features
It has 3 fields

4.1.1 속성을 통한 선택

GWR에서는 결측치가 있는 자료를 사용하지 못하므로 결측치가 없는 서울특별시 지역만 추출하였습니다.

# sel by attributes
`2017shp`@data$SIG_CD <- as.numeric(levels(`2017shp`@data$SIG_CD))[`2017shp`@data$SIG_CD]

seoul_shp <- subset(`2017shp`, SIG_CD <= 11740)

sigcd <- seoul_shp@data

4.1.2 인덱스 열 생성

위에서 분석할 때 사용하였던 데이터에는 시군구 코드가 없어서 조인이 제한되었습니다. 따라서 조인시 기준열을 만들기 위해 행정구역명의 형식을 같게 만들어주는 작업을 진행하였습니다.

a <- vector(mode = "character", length(sigcd$SIG_CD))

for (i in 1:length(sigcd$SIG_CD)) {
    if (as.numeric(sigcd$SIG_CD)[i] < 26000) {
        a[i] <- "서울특별시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 27000) {
        a[i] <- "부산광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 28000) {
        a[i] <- "대구광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 29000) {
        a[i] <- "인천광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 30000) {
        a[i] <- "광주광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 31000) {
        a[i] <- "대전광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 32000) {
        a[i] <- "울산광역시"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 37000) {
        a[i] <- ""  #세종
    } else if (as.numeric(sigcd$SIG_CD)[i] < 42000) {
        a[i] <- "경기도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 43000) {
        a[i] <- "강원도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 44000) {
        a[i] <- "충청북도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 45000) {
        a[i] <- "충청남도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 46000) {
        a[i] <- "전라북도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 47000) {
        a[i] <- "전라남도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 48000) {
        a[i] <- "경상북도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 49000) {
        a[i] <- "경상남도"
    } else if (as.numeric(sigcd$SIG_CD)[i] < 51000) {
        a[i] <- "제주특별자치도"
    } else {
        a[i] <- "NULL"
    }
}

do <- cbind(sigcd, a)

지역명 <- vector(mode = "character", length(do$SIG_CD))
for (m in 1:length(do$SIG_CD)) {
    지역명[m] <- paste(do$a[m], do$SIG_KOR_NM[m])
}
cbdo <- cbind(do, 지역명)

cbdo$지역명 <- as.character(cbdo$지역명)

finaldata <- inner_join(finaldata, cbdo)

map <- merge(seoul_shp, cbdo)

4.1.3 Join

공간 데이터프레임에 위에서 분석시 사용하였던 데이터를 조인하였습니다.

mapndata <- merge(map, finaldata)

4.1.4 셰이프 파일 도식

생성한 셰이프파일을 tmap을 사용하여 도식하였습니다.

tmap_mode("view")
tm_shape(mapndata) + tm_polygons("log.crime", id = "SIG_KOR_NM")

4.2 GWR

4.2.1 공간가중행렬 생성

추후 공간자기상관을 검정하기 위하여 공간가중행렬을 생성하였습니다.

xy <- coordinates(mapndata)

w <- poly2nb(mapndata, snap = 1)

class(w)

[1] "nb"

summary(w)

Neighbour list object:
Number of regions: 25 
Number of nonzero links: 110 
Percentage nonzero weights: 17.6 
Average number of links: 4.4 
Link number distribution:

2 3 4 5 6 7 
3 5 4 7 4 2 
3 least connected regions:
9 17 24 with 2 links
2 most connected regions:
2 18 with 7 links

plot(mapndata, col = "gray", border = "blue", lwd = 2)
plot(w, xy, col = "red", lwd = 2, add = TRUE)

wm <- nb2mat(w, style = "B")

datatable(wm, caption = "공간가중행렬", filter = "none")

4.2.2 커널 선정 및 GWR분석

최적 대역폯 산정 : Kernel 함수의 매개변수 최적값 산정방법은 교차검정법의 Score 함수법과 최소 AIC법으로 구분됩니다.

교차검증법

### CV method
gwr.1 <- gwr.sel(mapndata@data$log.crime ~ mapndata@data$log.GRDP + mapndata@data$police + 
    mapndata@data$log.cctv + mapndata@data$t.rate, mapndata)

Bandwidth: 14066.6 CV score: 0.7551871 
Bandwidth: 22737.51 CV score: 0.8064066 
Bandwidth: 8707.68 CV score: 0.7101516 
Bandwidth: 5395.686 CV score: 1.057663 
Bandwidth: 10754.6 CV score: 0.7208286 
Bandwidth: 7021.375 CV score: 0.7605218 
Bandwidth: 9489.535 CV score: 0.710756 
Bandwidth: 8988.553 CV score: 0.7095321 
Bandwidth: 9033.601 CV score: 0.7095331 
Bandwidth: 9009.305 CV score: 0.7095294 
Bandwidth: 9009.116 CV score: 0.7095294 
Bandwidth: 9009.035 CV score: 0.7095294 
Bandwidth: 9009.033 CV score: 0.7095294 
Bandwidth: 9009.033 CV score: 0.7095294 
Bandwidth: 9009.034 CV score: 0.7095294 
Bandwidth: 9009.033 CV score: 0.7095294

gwr.fit1 <- gwr(mapndata@data$log.crime ~ mapndata@data$log.GRDP + mapndata@data$police + 
    mapndata@data$log.cctv + mapndata@data$t.rate, data = mapndata, bandwidth = gwr.1, 
    se.fit = T, hatmatrix = T)

gwr.fit1

Call:
gwr(formula = mapndata@data$log.crime ~ mapndata@data$log.GRDP + 
    mapndata@data$police + mapndata@data$log.cctv + mapndata@data$t.rate, 
    data = mapndata, bandwidth = gwr.1, hatmatrix = T, se.fit = T)
Kernel function: gwr.Gauss 
Fixed bandwidth: 9009.033 
Summary of GWR coefficient estimates at data points:
                             Min.    1st Qu.     Median    3rd Qu.       Max.
X.Intercept.            2.6519751  3.6765541  4.6523765  5.2267136  6.7134392
mapndata.data.log.GRDP  0.0311756  0.0765880  0.1181178  0.1450739  0.2030228
mapndata.data.police   -0.0030252 -0.0026925 -0.0024245 -0.0021165 -0.0017549
mapndata.data.log.cctv -0.3523790 -0.2508988 -0.1913723 -0.1024642  0.0017622
mapndata.data.t.rate    2.3909962  6.3657456  6.8971908  7.9632365 13.8874641
                        Global
X.Intercept.            4.6045
mapndata.data.log.GRDP  0.1012
mapndata.data.police   -0.0023
mapndata.data.log.cctv -0.1234
mapndata.data.t.rate    6.1558
Number of data points: 25 
Effective number of parameters (residual: 2traceS - traceS'S): 11.78445 
Effective degrees of freedom (residual: 2traceS - traceS'S): 13.21555 
Sigma (residual: 2traceS - traceS'S): 0.1530389 
Effective number of parameters (model: traceS): 9.562252 
Effective degrees of freedom (model: traceS): 15.43775 
Sigma (model: traceS): 0.1415965 
Sigma (ML): 0.1112691 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 0.4574182 
AIC (GWR p. 96, eq. 4.22): -29.281 
Residual sum of squares: 0.3095203 
Quasi-global R2: 0.9168209

AIC

### AIC method
gwr.2 <- gwr.sel(mapndata@data$log.crime ~ mapndata@data$log.GRDP + mapndata@data$police + 
    mapndata@data$log.cctv + mapndata@data$t.rate, mapndata, method = "aic")

Bandwidth: 14066.6 AIC: -4.854602 
Bandwidth: 22737.51 AIC: -6.105851 
Bandwidth: 28096.43 AIC: -6.277339 
Bandwidth: 27415.83 AIC: -6.262654 
Bandwidth: 31408.42 AIC: -6.332105 
Bandwidth: 33455.35 AIC: -6.356084 
Bandwidth: 34720.42 AIC: -6.368288 
Bandwidth: 35502.27 AIC: -6.375019 
Bandwidth: 35985.49 AIC: -6.378905 
Bandwidth: 36284.13 AIC: -6.381209 
Bandwidth: 36468.7 AIC: -6.382598 
Bandwidth: 36582.77 AIC: -6.383443 
Bandwidth: 36653.27 AIC: -6.38396 
Bandwidth: 36696.84 AIC: -6.384278 
Bandwidth: 36723.77 AIC: -6.384474 
Bandwidth: 36740.41 AIC: -6.384595 
Bandwidth: 36750.7 AIC: -6.384669 
Bandwidth: 36757.06 AIC: -6.384715 
Bandwidth: 36760.98 AIC: -6.384744 
Bandwidth: 36763.41 AIC: -6.384761 
Bandwidth: 36764.91 AIC: -6.384772 
Bandwidth: 36765.84 AIC: -6.384779 
Bandwidth: 36766.41 AIC: -6.384783 
Bandwidth: 36766.77 AIC: -6.384785 
Bandwidth: 36766.99 AIC: -6.384787 
Bandwidth: 36767.12 AIC: -6.384788 
Bandwidth: 36767.21 AIC: -6.384789 
Bandwidth: 36767.26 AIC: -6.384789 
Bandwidth: 36767.29 AIC: -6.384789 
Bandwidth: 36767.31 AIC: -6.384789 
Bandwidth: 36767.32 AIC: -6.384789 
Bandwidth: 36767.33 AIC: -6.38479 
Bandwidth: 36767.33 AIC: -6.38479 
Bandwidth: 36767.34 AIC: -6.38479 
Bandwidth: 36767.34 AIC: -6.38479 
Bandwidth: 36767.34 AIC: -6.38479 
Bandwidth: 36767.34 AIC: -6.38479 
Bandwidth: 36767.34 AIC: -6.38479

gwr.fit2 <- gwr(mapndata@data$log.crime ~ mapndata@data$log.GRDP + mapndata@data$police + 
    mapndata@data$log.cctv + mapndata@data$t.rate, data = mapndata, bandwidth = gwr.2, 
    se.fit = T, hatmatrix = T)

gwr.fit2

Call:
gwr(formula = mapndata@data$log.crime ~ mapndata@data$log.GRDP + 
    mapndata@data$police + mapndata@data$log.cctv + mapndata@data$t.rate, 
    data = mapndata, bandwidth = gwr.2, hatmatrix = T, se.fit = T)
Kernel function: gwr.Gauss 
Fixed bandwidth: 36767.34 
Summary of GWR coefficient estimates at data points:
                             Min.    1st Qu.     Median    3rd Qu.       Max.
X.Intercept.            4.4916556  4.5623958  4.6213770  4.6633149  4.7618046
mapndata.data.log.GRDP  0.0956036  0.0986331  0.1013627  0.1037176  0.1069706
mapndata.data.police   -0.0023319 -0.0023041 -0.0022884 -0.0022660 -0.0022390
mapndata.data.log.cctv -0.1402540 -0.1328066 -0.1280813 -0.1229877 -0.1143196
mapndata.data.t.rate    5.9656521  6.0831367  6.2217449  6.3205499  6.4705119
                        Global
X.Intercept.            4.6045
mapndata.data.log.GRDP  0.1012
mapndata.data.police   -0.0023
mapndata.data.log.cctv -0.1234
mapndata.data.t.rate    6.1558
Number of data points: 25 
Effective number of parameters (residual: 2traceS - traceS'S): 5.568992 
Effective degrees of freedom (residual: 2traceS - traceS'S): 19.43101 
Sigma (residual: 2traceS - traceS'S): 0.1693034 
Effective number of parameters (model: traceS): 5.293462 
Effective degrees of freedom (model: traceS): 19.70654 
Sigma (model: traceS): 0.1681156 
Sigma (ML): 0.14926 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): -6.38479 
AIC (GWR p. 96, eq. 4.22): -18.8629 
Residual sum of squares: 0.5569633 
Quasi-global R2: 0.8503241

Adaptive kernel

### adaptive gwr.adapt <- gwr.sel(mapndata@data$log.crime ~ mapndata@data$log.GRDP
### + mapndata@data$police + mapndata@data$log.cctv + mapndata@data$t.rate,
### mapndata, adapt = TRUE)

# gwr.fit3 <- gwr(mapndata@data$log.crime ~ mapndata@data$log.GRDP +
# mapndata@data$police + mapndata@data$log.cctv + mapndata@data$t.rate, data =
# mapndata, adapt = gwr.adapt, se.fit=T, hatmatrix=T)

4.2.3 GWR 모형 분석결과

분석결과, CV에 의한 Band Width \(\theta \ = 9009.033\)은 잔차제곱합 0.3095203, \(R^2 = 0.9168209\), AIC에 의한 Band Width \(\theta \ = 36767.34\)은 잔차제곱합 0.5569633, \(R^2 = 0.8503241\)로서, 교차검증법이 적합하다고 판단되어 gwr.fit1 모델을 채택하였습니다.

gwr.fit1

Call:
gwr(formula = mapndata@data$log.crime ~ mapndata@data$log.GRDP + 
    mapndata@data$police + mapndata@data$log.cctv + mapndata@data$t.rate, 
    data = mapndata, bandwidth = gwr.1, hatmatrix = T, se.fit = T)
Kernel function: gwr.Gauss 
Fixed bandwidth: 9009.033 
Summary of GWR coefficient estimates at data points:
                             Min.    1st Qu.     Median    3rd Qu.       Max.
X.Intercept.            2.6519751  3.6765541  4.6523765  5.2267136  6.7134392
mapndata.data.log.GRDP  0.0311756  0.0765880  0.1181178  0.1450739  0.2030228
mapndata.data.police   -0.0030252 -0.0026925 -0.0024245 -0.0021165 -0.0017549
mapndata.data.log.cctv -0.3523790 -0.2508988 -0.1913723 -0.1024642  0.0017622
mapndata.data.t.rate    2.3909962  6.3657456  6.8971908  7.9632365 13.8874641
                        Global
X.Intercept.            4.6045
mapndata.data.log.GRDP  0.1012
mapndata.data.police   -0.0023
mapndata.data.log.cctv -0.1234
mapndata.data.t.rate    6.1558
Number of data points: 25 
Effective number of parameters (residual: 2traceS - traceS'S): 11.78445 
Effective degrees of freedom (residual: 2traceS - traceS'S): 13.21555 
Sigma (residual: 2traceS - traceS'S): 0.1530389 
Effective number of parameters (model: traceS): 9.562252 
Effective degrees of freedom (model: traceS): 15.43775 
Sigma (model: traceS): 0.1415965 
Sigma (ML): 0.1112691 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 0.4574182 
AIC (GWR p. 96, eq. 4.22): -29.281 
Residual sum of squares: 0.3095203 
Quasi-global R2: 0.9168209

4.2.4 GWR 모형의 Moran Test

지리가중회귀 모형에 대한 Moran Test 결과는 다음과 같이 요약됩니다.

gwr.morantest(gwr.fit1, nb2listw(w))


    Leung et al. 2000 three moment approximation for Moran's I

data:  GWR residuals
statistic = 11.966, df = 15.663, p-value = 0.7255
sample estimates:
         I 
-0.2087221

검정통계량은 11.966, \(P=0.7255 > \alpha \ = 0.05\) 로서 공간적 자기상관이 없다는 귀무가설 기각을 실패합니다. 즉, 지리적 가중회귀 모형의 Moran’s I = -0.208로서 강한 공간적 자기상관은 없는 것으로 판단됩니다.

4.2.5 분석결과 plot

모형을 토대로 분석을 수행한 분석결과를 도시해보았습니다. 전역적 모델인 OLS모델보다 국지적 모델인 GWR의 잔차제곱합이 가장 작았습니다.

### pred vs actual

lmresd <- gwr.fit1$lm$residuals
lmfit <- gwr.fit1$lm$fitted.values

seoul_shp1 <- seoul_shp
seoul_shp1@data <- cbind(seoul_shp@data, lmresd, lmfit)

실제 값과 추정 값

지도

실제 값

tm_shape(mapndata) + tm_polygons("log.crime", id = "log.crime", alpha = 0.6)

OLS

tm_shape(seoul_shp1) + tm_polygons("lmfit", id = "lmfit", alpha = 0.6)

GWR

tm_shape(gwr.fit1$SDF) + tm_polygons("pred", id = "pred", alpha = 0.6)

표

value <- cbind(as.character(mapndata@data$SIG_KOR_NM), mapndata@data$log.crime, seoul_shp1@data$lmfit, 
    gwr.fit1$SDF$pred) %>% as.data.frame()

names(value) <- c("구", "Actual", "OLS", "GWR")

datatable(value, caption = "Actual vs Pred", filter = "top")

OLS와 GWR의 잔차 비교

지도

OLS

### 잔차
tm_shape(gwr.fit1$SDF) + tm_polygons("gwr.e", id = "gwr.e", alpha = 0.6)

GWR

tm_shape(seoul_shp1) + tm_polygons("lmresd", id = "lmresd", alpha = 0.6)

잔차 표

value1 <- cbind(as.character(mapndata@data$SIG_KOR_NM), seoul_shp1@data$lmresd, gwr.fit1$SDF$gwr.e) %>% 
    as.data.frame()

names(value1) <- c("구", "OLS", "GWR")

datatable(value1, caption = "OLS vs GWR residuals", filter = "top")

SSE 표

osse <- vector(mode = "numeric", 25)

OLS <- as.numeric(levels(value1$OLS))[value1$OLS]

for (m in 1:25) {
    osse[m] <- OLS[m] + 1
}

OLS_SSE <- sum(osse)

gsse <- vector(mode = "numeric", 25)

GWR <- as.numeric(levels(value1$GWR))[value1$GWR]

for (n in 1:25) {
    gsse[n] <- (GWR[n])^2
}

GWR_SSE <- sum(gsse)

sse <- cbind(OLS_SSE, GWR_SSE)

sse

     OLS_SSE   GWR_SSE
[1,]      25 0.3095203

GWR Local R2값 확인

### localr2
tm_shape(gwr.fit1$SDF) + tm_polygons("localR2", id = "localR2", alpha = 0.6)

5 Conclusion

범죄율은 다양한 사회적, 공간적 특징에 의해 영향을 주고받는 요인으로서 여러 지리학의 분야에서 연구대상이 되어오고 있습니다. 저희 조는 범죄율에 영향을 미칠 수 있는 요소들을 선택해 여러 통계적 분석을 진행하였습니다.
먼저, 인구 1만명당 CCTV와 경찰의 수, GRDP, 순이동, 그리고 기초 수급자의 비율을 범죄율에 영향을 미칠 것으로 예상되는 요소들로 선정하였습니다. 이 데이터들을 정제한 뒤, 특정 요소에 대한 단순선형회귀 분석을 하였습니다.
그 결과, GRDP와 범죄율의 로그값은 강한 양의 상관관계를 보였으며 이를 유의미한 분석으로 확인하였습니다. 다음으로는 다중회귀분석을 집중적으로 진행하였고, 전진선택법과 후진제거법 등의 stepwise regression을 활용해 분석한 결과 GRDP, CCTV의 수, 기초수급자의 비율은 범죄율과 양의 상관관계를, 순이동과 치안력은 음의 상관관계를 보였습니다. 이 중 GRDP, 치안력, 기초수급자의 비율이 범죄율에 상대적으로 큰 영향을 미치는 요소들로 드러났습니다. 이후 잔차분석을 실행해 독립성, 선형성, 다중공선성을 검정하였고 분석에 이상이 없음을 확인했습니다.
이어서 Influence plots와 outlier test 등에서 outlier로 판명된 103번을 제거하여 다시 회귀분석을 진행하였고, 이로 인해 변수 간 상관관계가 일부 강화되었습니다. 추가적으로 지리가중회귀분석을 진행한 결과 선형회귀분석에 비하여 잔차가 줄어든 것을 확인 할 수 있었습니다.

공간정보분석

[ASS-06] Team Zulu ‘Attack Launched!’ by: 서용훈, 김예원, 김예진, 박경준